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DEFINITION OF SYMBOLS 


Standard indicial and matrix notations are used throughout this paper. 
Repeated indices, unless enclosed by parenthesis, indicate summation. Upper- 
case Latin indices generally indicate points in space, whereas lower-case Latin 
indices indicate elements of an array. Greek indices are, in general, associated 
with the in-plane coordinate system and range from i to 2. The following symbols 
are used. 
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exp 


„N 
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-m 
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NO! 
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a 


P Ni 


Determinant of a _ 

ap 

Coefficients of the first fundamental form of the undeformed 
surface 

Tangent base vector in the undeformed middle plane 

Element identification number 

Bilinear nodal function 

Metric tensor in the undeformed state 

Thickness of plate 

Bending stiffness matrix 

Membrane stiffness matrix 

Shear stiffness matrix 

Total number of assembled nodes in system 

Local nodal moments 

Unit vector normal to middle plane 

Surface force components 

Local nodal loads 
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DEFINITION OF SYMBOLS (Continued) 
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N 

U Surface displacement vector 

U General displacement vector 
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X^ T Surface coordinates of node N 

No 

Z Normal to the surface coordinate 
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ij 
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Na 

P Point on middle surface of undeformed plate 
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P Body force 

P^ T . Global nodal forces 

Ni 

U Total strain energy 

U q Strain energy density 
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Global in-plane displacements 

General vector field 
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Global value of V at node N 

Global transverse displacements 
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A STUDY OF STIFFNESS MATRICES FOR THE 
ANALYSIS OF FLAT PLATES 

SUMMARY 


The analysis of thin plates in bending is considered with four different 
rectangular finite element representations. The first representation approxi- 
mates the transverse displacement by a sixth-order two-dimensional generaliza- 
tion of a Hermitian interpolation polynomial. This representation requires that 
the Kirchhoff hypothesis be satisfied throughout the element. Therefore the 
transverse shear strains vanish, and the element exhibits no shear stiffness. 

The second representation uses a simple bilinear approximation for the trans- 
verse displacements and rotations without the Kirchhoff hypothesis. In the latter 
case both shear-stiffness and bending-stiffness matrices are obtained. The third 
approximation uses a discrete Kirchhoff hypothesis which introduces a constraint 
between the transverse nodal displacements and nodal rotations. This case also 
contains both bending and shear stiffness. The fourth representation uses only 
the bending stiffness developed in the previous case. This is a logical approxi- 
mation since the discrete Kirchhoff hypothesis causes the transverse shears to 
vanish in the limit. Numerical examples are presented to demonstrate the 
relative accuracy of the finite elements investigated. 


INTRODUCTION 
General Comments 

Since the thickness of a plate is small in comparison with other dimen- 
sions , certain simplifying assumptions can be introduced which reduce plate 
problems to two-dimensional rather than three-dimensional analysis. The well- 
known Kirchhoff hypothesis is an example; it assumes that lines normal to a 
plate's middle surface before deformation remain normal after deformation. 



Another less restrictive assumption asserts that displacements vary linearly 
over the plate thickness. All of these assumptions are kinematic in nature; they 
impose no restrictions on the order of magnitude of the strains in the plane of 
the plate or on the order of magnitude of the displacements. 

When a thin flexible plate is subjected to transverse loads, it displaces 
normal to its middle plane and forms a curved surface. If the transverse dis- 
placements are small in comparison with the thickness of the plate , the strains 
in the middle surface are usually small and negligible in comparison with those 
developed in the extreme fibers. If the plate undergoes large transverse dis- 
placements, however, significant strains may be developed in the deformed 
middle plane. On the other hand, if a flexible plate is subjected to sufficiently 
large loading in its plane, it will buckle laterally, and bending stresses will be 
developed. Mathematical descriptions of these phenomena involve highly non- 
linear partial differential equations in the transverse displacements. Few exact 
solutions to these equations are available in the literature, and, in the case of 
plates with irregular shapes and boundary conditions, exact solutions are prac- 
tically intractable even when classical linear theory is used. Because of this 
difficulty in obtaining exact solutions , numerical methods are often employed to 
obtain quantitative solutions to these problems. 

Among the numerical methods available , the finite element method is 
appealing and is the technique investigated in this thesis. In the finite element 
method, a continuous plate is represented by an assembly of small polyhedral 
plate elements of finite dimension, each of which is assumed to have finite de- 
grees of freedom. The displacement fields within each element are approxi- 
mated by polynomial functions of the local coordinate system associated with 
each finite element. These discrete structural elements are interconnected at 
a finite number of node points. The problem is then reduced to one of determin- 
ing a finite number of unknowns. 


Scope of Paper 


Ironically the finite element formulation of plate problems involves some 
complications not encountered in finite element analyses of three-dimensional 
bodies. In the three-dimensional body formulation, only the relative displace- 
ment of each node point is of concern, whereas in the thin flat-plate formulation, 
in addition to the relative displacement of the nodes , the relative rotations and 
twist come into play. Because of these complications most finite element for- 
mulations for the analysis of plates involve high-order polynomial approximations 
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for the transverse displacement field. These approximations prove unwieldy 
when extending the formulations to shell analyses and to geometrically nonlinear 
plate-problems. 

Several investigations [1-9] have developed linear finite element stiffness 
matrices for the analysis of thin plates in bending. Of particular interest among 
these is the paper by Clough and Tocher [2] which investigates the relative ac- 
curacy of seven different types of finite element representations. Melosh*s paper 
[4] was the first to utilize a rectangular- shaped finite element representation 
for the analysis. 

Bogner, Fox, and Schmidt [1] presented interpolation formulas in 
orthogonal curvilinear coordinates, which they used to approximate the displace- 
ment field within a finite element of a flat rectangular plate. These interpolation 
formulas are polynomials of sixth order in the coordinates. Although this re- 
presentation is highly complex, it yields a stiffness matrix which exhibits good 
convergence characteristics. This matrix is examined later in this paper. 

Utku [ 8] and Melosh and Utku [ 9] have developed a triangular discrete 
element formulation adaptable to shells. This representation utilized simple 
linear approximations for the displacement and rotation fields. In order to 
obtain convergence, however, it was necessary to modify the stiffness matrix 
acquired in this formulation. Several different schemes were utilized to modify 
the matrix. This approach was abandoned in the present study in hope of obtain- 
ing a more rational approach to the problem. 

It appears, therefore, that finite element solutions to plate problems in 
the past have been extremely involved because of the Kirchhoff hypothesis and 
convergence criteria involving continuity requirements of the slopes. A more 
simple finite element representation is needed before the method can be extended 
to general shell problems and nonlinear plate problems. 

This paper, therefore, investigates several finite elements based on 
relatively simple displacement approximations which satisfy the convergence 
criteria, but with one exception do not satisfy the Kirchhoff hypothesis through- 
out the element. The purpose of the study is to develop a simple finite element 
model which converges well enough. The discrete variables selected in this 
study are the displacements and rotations of the nodal points. Although this 
paper considers only the classical small deflection theory of plates, the finite 
element representations presented are ones which lend themself readily to the 
shell analysis and to nonlinear plate problems. Numerical examples are pre- 
sented to demonstrate the relative accuracy of the finite elements investigated. 
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KINEMATIC CONSIDERATIONS 
The Discrete Model 


In finite element formulations of plate problems, a continuous plate, such 
as that shown in Figure la, is represented by an assembly of a finite number E 
of small plate elements as indicated in Figure ib. The geometry of a typical 
element is defined by the location of a number of nodal points on the element's 
boundaries, connected by smooth curves called nodal lines. In contrast to 
classical analyses wherein relationships between mean values of certain varia- 
bles associated with differential elements are obtained and the dimensions of the 
elements are shrunk to zero as their number becomes infinite, the dimensions of 
finite elements remain finite throughout the analysis. Ideally these dimensions 
are small in comparison with characteristic dimensions of the assembled system 
so that displacement fields within each element are adequately approximated by 
appropriate functions (usually polynomials) of the coordinates. 




FIGURE 1. FINITE ELEMENT 
REPRESENTATION 


In this study the locations of 
points in the discrete system are given 
by the curvilinear surface coordinate 

Q> 

X (a a 1, 2) which are embedded in 
the middle surface of the plate and a 
coordinate Z normal to the middle sur- 
face. The undeformed shape of the 
middle surface of each finite element 
is assumed to be a flat, curvilinear 
quadrilateral, and for convenience the 
vertices of these quadrilaterals are 
selected as the nodal points of the ele- 
ments ( Fig. ib) . 

A vector field V(X° ! , Z) defined 
throughout the domain D occupied by the 
distributed plate, identifies with every 
point P in D a vector V( P) . In the finite 
element representation of D, the field 
V is depicted by a finite set of quantities 
which represent the values of V at each 
nodal point, the values of V at other 
points being given by appropriate inter- 
polation formulas. Thus, if there are 
m nodal points in the discrete model 
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after all elements have been connected together to form a single assembled 
system and if denotes the value of V at node N, then the set V^(N = i, 2, 

. . . m) is the global representation of the field V. On the other hand, if of the 

totality of E finite elements element e is isolated and examined independently 

of the other elements , then the field V is characteristized within e by the set 

(M = 1, 2, 3, 4) , where V,, is the value of V at node M of element e. The 
~Me —Me ~ 

sets of quantities = 2 ’ 3, 4; e = 1, 2, 3, ... E) are referred to as 

local representations of V corresponding to elements e = 1, 2, 3, ...E. 

This distinction between global values V XT and local values V,, of the 

same field V is introduced for convenience. The general behavior of a typical 
finite element can be described in terms of local values of various fields, in- 
dependent of its mode of connection with or the behavior of adjacent elements. 
The local and global values are then related through transformations of the form 


V = £2 V 
-Me MNe ~N 


( 1 ) 


wherein M = 1, 2, 3, 4; N = 1, 2 . . . m; e = 1, 2 . . . E 


! 1 if node M of element e is identical to node N of the 

assembled system ( 2) 

0 if otherwise. 

The transformation defined in equation ( 1) is said to establish the connectivity 
of the discrete system. Mathematically it establishes the required dependencies 
between local and global values of the field V; physically it connects the elements 
together at their nodes to form a single unit. 


Bilinear Vector Field Approximation 


This discussion is limited to a typical element of the system. For 
convenience the element index e is temporarily dropped. 

Assume that the components of a vector field within the region of the 
rectangular finite element can be approximated by the form 
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( 3 ) 


y a = B a + C a + D_ X*X 2 

r * 


a 


01 01 01 

where B , C ^ , and D are the undetermined constants. 


oi oi 

Let V N and X^. ( N =1,2, 3, 4; a = 1, 2) denote respectively the com- 
ponents of the vector field at node N and the surface coordinates of node N of a 
typical finite element. To obtain the undetermined constants in terms of nodal 
quantities, equation ( 3) is evaluated at each of the four node points of the ele- 
ment being considered. Evaluating equation ( 3) for the component V 1 


V 1 = 

N 


B 1 + C 1 


x * + 

P X N 


D 1 Y 


N 


(4) 


where 


y n = x1 (n) x2 (n) ( nosumonN) 


In matrix form 




Cb 


where 


= {v\, v|, vl, vi) 

b = {B 1 , Cl, Cl, D 1 } 


and 

1 x{ x\ Y 1 
1 x| x| y 2 
1 x! x| y 3 
1 xl x 2 4 y 4 


(5) 


( 6 ) 


(7a) 

(7b) 


( 8 ) 
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Solving equation ( 6) for b , 


b = C' 1 vA. 

'N 




(9) 


Treating the other component of the vector field in a similar fashion, we 
find that 

OL N ca 

b = k v“ 

N 


N _.Q! 

°f> = °f> V N 

OL N CL 

D = d V 

N 


( 10 ) 


where 




N pNMRS \ 

°n “xe s “m x e y s 


(ID 


and 




„NMRS _ J_ 

6 C e NMRS 


( 12 ) 


In these equations I s the four-dimensional permutation symbol, C is 

the determinant of C, and a^ = 1(M = 1, 2, 3, 4) . 

Substituting equation (10) into equation ( 3) gives 


a 


N 


r a 

'N 


N 

“/3 


OL 

7 

N 


V" = k*' V„ T + cl’ X'' vZ + d N X 1 X 2 V 


a 

N 


(13) 


Further, let 

f N .k N + c*x fS + d N x‘x=, 


then 


(14) 
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( 15 ) 


a N a 
V = f V 

N 


Equation ( 15) gives the bilinear approximation for the vector field component 

OL N 

V in terms of the nodal values of this vector field. Note that f is independent 

OL 

of the direction of the vector component V . 


Lagrangian Strain Tensor 


In order to obtain a finite element representation for a thin plate , it is 
convenient to review briefly the kinematics of thin shells [ 10] . In the following, 
Latin indices range from 1 to 3. 

The general definition of the Lagrangian strain tensor y .. in the case of a 
three-dimensional continuum is ” 


y.. = "5“ (G.- ~ g..) 

i] 2 ij to ij 


( 16 ) 


where G.. and g.. are the metric tensors in the deformed and undeformed states 
ij 

respectively. If G is expressed in terms of g and derivatives of the displace- 

men field U and if linear small deflection theory is used, then equation ( 16) 
reduces to the linear strain-displacement relations 


y.. = ~(u. . + u. .) 

IJ 2 i;j j;i 


( 17) 


where U are the covariant components of the displacement vector and the 
semicolon denotes covariant differentiation with respect to a curvilinear system 

OL 

of convected coordinates X (a = 1, 2). 


Several simplifications of equations ( 16) and ( 17) are possible in the 
case of finite deformations of continuous flat plates. Consider, for example, 
the initially flat thin plate shown in Figure 2, the geometry of which is described 

OL 

by the curvilinear surface coordinates X ( a - i , 2) and a coordinate Z normal to 
the middle surface. The position vector of a general point P in the undeformed 
plate is denoted r. Tangent base vectors in the undeformed middle plane are 
denoted a and n denotes a unit normal to the middle plane. 

~Q! ~ 
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X' 


FIGURE 2. GEOMETRY OF DEFORMATION 

If the plate undergoes a general deformation, points P and P move with 
displacement vectors U and U to locations P* and P* . In the deformed state 

— rv rv/ — 

point P* is located by the position vector R, and point P* in the deformed 
middle surface is located by the position vector R, as shown in Figure 2. The 
vector initially normal to the middle surface is rotated after deformation. This 
rotation vector is denoted 0. From Figure 2 it is clearly seen that 

rv/ 

U = U + Z(n x 0) (18) 

rj r%j rv rv/ 


and 
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(19) 


R = R + Z N = r + U + Z(n + n x ©) 


rsj r^J rsj r+-> r+j 


r+j n*J 


where N denotes a unit vector tangent to the deformed Z coordinate line. The 
metric tensor in the deformed state is given by 


G 0 — R, • R, „ . 

Oi j3 ~ a ~ j8 


( 20 ) 


Combining equations (20) and (19) and neglecting nonlinear terms in U and © 
as well as a term involving Z 2 , it is found that 


^ a p a a (3 + p + Z *a/3 


( 21 ) 


where 


a . = a • a_ 
ap ~oi ~(3 


2y a» ‘ U a;/S + U /3;c 


( 22 ) 

(23) 


and 


X R = 's/a" ( e R Q X + e © A ) . 
aft XP ;a \a ;P 


(24) 


In the above equation the determinant of a^ is denoted as a. Similarly, 


G = R, • N = w, + */a" e © 
<23 ~ <2 ~ a <2A 


(25) 


Finally, substituting equations (21) through (25) into equations (16) yields 


V = T <*0 + Z V 


(26a) 


and 


y = w, + <s/a" e , © A . 
r a3 ex a\ 


( 26b) 


Note that the transverse extensional strain y g3 was assumed negligible 
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Kinematics of the Discrete System 


To complete the kinematic formulation, it is necessary to obtain the 
strain displacement relations developed in the preceding section in terms of 
nodal displacements and rotations. The components of in-plane displacements 
and rotations may be approximated by a bilinear form. A higher-order form for 
w, however, will also be investigated. Thus 


N 

U = f U T 
a Na 


(27a) 


e a -f"e 


Na 


(27b) 


N 

but w = w \v T 
T N 

N N 

where the special case of ij) = f is one case to be investigated. 


(27c) 


and 


Combining equations (27) with equations (23) , (24) , and (26) , 

y = — (f N U + f N U \ + Z\Ta (e f N ©* + e f N 0* 
'a/3 2 \ \a N/3 ;j3 N a I \ 7^3 ;a N \cx ;[3 b i 


N i— £ N q A 

y = lb W + N/a € , f 0^ T 
r ce3 r , a N a X N 


) 


(28a) 

(28b) 


These equations express the strain at any point in a typical plate element 
in terms of the nodal displacements and nodal rotations. It has been pointed out 
that the global and local values of a vector field could be related through trans- 
formations of the type given by equation ( i) . Applying these connectivity trans- 
formations to the nodal displacements and rotations. 


U = S2 U 

Mae MNe Na 

(29a) 

W Me = °MNe W N 

(29b) 

e“ e 01 

Me MNe N 

(29c) 
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where U^, and ©^ T are the global values of the in-plane displacements, 

transverse displacement and rotations at global node N, respectively. The con- 
nectivity matrix was defined previously by equation (2) . 

Combining equations (28) and (29) , 


-IN N 

y - ± ( f* ft U + f O U ) 

r a(3e 2 ; 0 i NMe M/3 ;/3 Me Ma' 


y = i/) N S2 W + 's/a" e f ^ £2 0 ^ . 

y ot 3e v ,ct NMe M a\ NMe M 


„N 


(30a) 

(30b) 


Equations ( 30) express the strain in a typical element e of the discrete 
system in terms of the global nodal displacements and rotations. 


STIFFNESS RELATIONS 
Local Stiffness Relations 


It is possible to develop the stiffness relations for a typical finite element 
of a thin plate from energy considerations. The total potential energy in an 
elastic body is given by 

n = U + 0 (31) 

where U = J U dV is the total strain energy of the finite element and is the 
v 

potential energy of the extrinsic forces acting on the body. The potential energy 
of the extrinsic forces may be written as 

S2 = - f P wdV - f p a U dS t (32) 

J J 01 

V Sj 

OL 

where p is the prescribed body force acting in the Z direction and p are the 
components of the prescribed surface force vectors. 
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The integration of the surface forces is taken over the area Sj, on which 
the forces are prescribed. From equation ( 18) it follows that 


U = U + Zn/T e . 
a a a/3 


(33) 


The energy function £2 may now be written in terms of nodal displacements and 
rotations by combining equations ( 32) , ( 33) , and the approximation equations 
( 27) . 


« = - /p sla dVw N - / p a f N dSj U Na 

V Sj 

- / p“ Z^r £ f dS, e£ (34) 

Si 

The nodal quantities have been taken outside the integrals since they are 

a 

independent of the coordinates X , Z. The integral terms in this equation are 
components of the so-called consistent load vectors and are denoted as 


Na r a N __ 

p = J p f dS t 

(35a) 

si 


N3 r N , — 

p = J p ip nT a dV 

(35b) 

V 


N r Of i — 

m = I p Z\[a e . f dSn . 

a/3 1 

(35c) 


si 


Ni 

where p (i = i,2,3;N=i,2,3,4) represents the generalized nodal force 

N 

acting at node N of a typical element and m (a = 1,2) represents the generalized 
nodal moment at node N. 

Assuming that the deformation is reversible, either isothermal or 
adiabatic, an elastic potential function exists which represents the strain 

energy per unit volume of the undeformed element. The strain energy function 
can be written in terms of the strains. By means of equations ( 28) , it can also 
be written in terms of the nodal displacements and rotations. The strain energy 
density function for an elastic flat plate element is of the form 
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( 36 ) 


TT 1 ijkl - - 

U o = 2" E y ij y kl 


4jkl 


where E J are elastic constants. The total potential energy is then 

n = f/ E ‘ ! “ Vki dv - pN “ u 


Na 


(37) 


N3 N 

-p w XT - m 9 . 

N a N 

According to the principle of minimum potential energy the strained 
element reaches an equilibrium state when 


9ii _ an = an 

au a 

Na 8w xt ae XT 
N N 


(38) 


Thus from equation (37) it follows that 


Na 

P 


ay. 


/ E 


ijkl- 


kl 


%. au XT 

ij Na 


dV 


(39a) 


N3 


= f E 


ijkl - 8y kl 



ij 9w 


dV 


N 


(39b) 


m 


N 

' a 


= f E 


ijkl - 8y kl 
y ij ae« 


dV 


(39c) 


Since y. are linear functions of the nodal displacements and rotations, the right 
side of equations (39) are linear functions of these nodal quantities. 

Equations ( 39) may be rewritten in more convenient matrix notation. 


P = k 6 (40) 

r>~> r^> * 

where 
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2 * (p 11 , p 21 > • • • p 41 > p 12 . • • • p 43 ; m l. “a 2 . • • • m l} 


(41a) 


and 


6 = {Un, U 21 , . . . U 41 ; U 12 , ... U c ; w t , ... w 4 ; ©|, ... ©|; © 2 4 . . . © 2 4 } . 

(4ib) 


k is a symmetric matrix which expresses the linear relationship between 
the generalized nodal forces and displacements for a typical element. It is 
convenient to express k by the sum 


k = k + k. + k 
~ ~m ~D ~ s 


(42) 


where 


k = 
m 


r ijkl - 9 ^kl 

J E J y.. ^ — dV 


ij au 


(43a) 


nee 


, r x-.ykl - ^kl , TI 

k, = I E J y.. dV 

~h J n 


ij 8w 


(43b) 


N 


k , j E nn - 

~S J 1] 

v J 80 

N 


dV 


(43c) 


The k , k, , and k matrices are known as membrane, bending, and shear 
~m ~b ~s 

stiffness matrices, respectively. 


Global Stiffness Relations 


The stiffness relations derived in the previous section describe the 
elastic behavior of a single typical rectangular flat plate element relative to a 
local reference frame. They are independent of the location of the element in 
the assembled system, the boundary conditions, and the loading conditions. To 
describe the behavior of a given structure with a specific shape and specific 
boundary conditions, it is necessary to assemble the elements. 
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It is convenient to select the local coordinate systems associated with 
each rectangular element of the assembled system parallel to the global coordi- 
nate system. Then it is not necessary to rotate the local nodal vector quantities 
to the global reference frame. 


can be 
that 


The transformations of the previous Section, Kinematic Considerations, 

conveniently rewritten in matrix notation by introducing a matrix £2 such 

^ 0 


6 = £2 A (44a) 

~e ^6 r '“' 


and 


P = £2 p 


(44b) 


where 


A = {Uu, TJ a , Wj, ©i, 0 2 i; ... U U w , e 1 ,© 2 } 
~ u 11 u 1 1 1 ml m2 m m m J 


(45a) 


and 


P = {P U , P 12 , P 13 , M}, M 2 1 , 


ml m2 m3 , m „ m, . A , . 

P , P ,P , Mj , Mj} . (45b) 


P and are the global values of generalized forces and moments 

at global node N, respectively. External generalized node forces and generalized 
node displacements of the assembled system are related by the global stiffness 
relation 


P = K A . (46) 

Combining equations (40) , (44) , and (46) , 

E T 

K= V £2 k £2 . (47) 

r*J L — > 0 < — ' 0 r^0 

e = l 


The matrix K is known as the total unsupported stiffness matrix for the assembled 
system. Note that the connectivity relations in equation (47) do not take into 
account any rotation of coordinate systems since it is assumed that the local and 
global coordinate systems are parallel. 
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Boundary conditions are next applied to the assembled system by pre- 
scribing generalized displacements or forces at appropriate nodes. In this way 
the final supported stiffness matrix is obtained from the unsupported matrix. 
Then equation ( 46) reduces to a system of independent linear algebraic equations 
in the unknown nodal displacements and rotations. Once these are solved, the 
strains are obtained for each element from equations ( 30) . The element 

stresses o’ 1 "* may be evaluated by applying the linear stress strain relation 

G 



e 


= E 



(48) 


This concludes the formulation of the problem. Several special cases are 
considered in the following section. 


PLATE ANALYSES 
General Comments 


All of the plate stiffness matrices considered in this discussion were 
developed using the procedure described in the previous Sections specialized to 
rectangular elements and Cartesian coordinates. Displacement and rotation 
fields are approximated by relatively simple functions which satisfy the conver- 
gence criteria but with one exception violate the Kirchhoff hypothesis of normals 
remaining normal throughout the element. The types of deformation and rotation 
approximations on which each of the analyses is based will be discussed briefly. 
The resulting bending, shear, and membrane-stiffness matrices are presented 
for each case. The coefficients of these stiffness matrices are listed in the 
Appendix. 


Special Cases 


As mentioned in the preceding Section, Stiffness Relations, the elemental 
stiffness matrix may be expressed as a sum of the membrane, bending, and 
shear matrices. 

k = k + k, + k . (49) 

~ ~m ~b ~s 
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The membrance stiffness matrix depends on the in-plane displacements 
but is independent of the transverse displacement w and rotations ©^ . On the 

other hand, the in-plane displacements do not influence the bending and shear 
stiffnesses. The in-plane middle surface displacements are approximated by a 
bilinear form 


U 

a 



U-, • 

No; 


(50) 


The strain-displacement relations for a typical element with membrance strains 
only is of the form 




+ U 


P,a 


(51) 


The membrane-stiffness relation is obtained by combining equations (51) and 
( 39a) and is of the form 


p = k U 
& m ~m ~m 


(52) 


where 


Pi 2 > P 2 i> P 22 > Psi » P 32 > P4l> P 42 } 


(53a) 


and 


U m = {Uii, U 12> U 21) U 22 , U 31 U 32 , U 41 , U 42 } . (53b) 

The coefficients of the membrane stiffness matrix are listed in Table I. 
All of the cases considered in this report assume the same bilinear form for U , 

and therefore they all have the same membrane stiffness matrix. This matrix 
has been presented frequently in the literature [1, 7] and is known to yield good 
results. Therefore, membrane stiffness matrices will not be discussed further. 

Case I. The vector field approximation to be utilized in the following 
cases is of the simple polynomial form. As mentioned earlier, the bilinear 
approximation is mathematically simpler than the higher order polynomial ap- 
proximations generally used; however, the relative accuracy obtained from the 
two different approximations should be compared. Toward this end a higher 
order polynomial approximation for the transverse displacement field is con- 
sidered in this case. 
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TABLE I. MEMBRANE STIFFNESS MATRIX 


Row 

Column 


1 2 

3 

4 

5 

6 

7 

8 

1 

2(bj + b 2 ) fc>3 

( b t - 2b 2 ) 


-<b,+ b;,) 

-bs 

( ~2b 1 + b 2 ) 

-hi 

2 

2(b 5 + b 6 ) 

-b 4 

(2b 5 + b 6 ) 


-(b 5 + b 6 ) 

b 4 

( b 5 - 2b e ) 

3 


2(b 1 + b 2 ) 

-bs 

( -2bj + b 2 ) 

b 4 

-<b, + b 2 ) 

b 3 ' 

4 



2(b 5 + b e ) 

-b 4 

(b 5 - 2bg) 

b 3 

-(b 5 + b 6 ) 

5 




2(b 4 + b 2 ) 

b 3 

(bj - 2b 2 ) 

b 4 

6 

SYMMETRIC 




2(b g + b 6 ) 

-b 4 

( 2b 5 + b e ) 

7 






2(b t + b 2 ) 

-b» 

8 







2( b 5 + b e ) 


Note: Symbols used in this table are defined in the Appendix. 


To simplify the following derivation, fibers initially normal to the middle 
surface are assumed to remain normal after deformation; that is, the Kirchhoff 
hypothesis is invoked. Hence the shear strain components are required to 
vanish. The displacement of any point in the plate is then given by 

U = U - Zw (54) 

a a ,a 


The strain displacement relations for this case become 
2y =U . + U_ - 2Zw 


aj3 a, j3 (3, a 


,a/3 


(55) 


Turning to the representation of w, it is important first to note that it is 
not sufficient merely to match values of w along the nodal lines to assure 
monotonic convergence [5]. Because of the Kirchhoff assumptions, plate de- 
formation is dependent on displacements, slopes, twists, curvatures, and so on 
of the middle surface. Consequently, the displacements of points not on the 
middle surface represent at best mean values subject to the initial requirement 
that normals remain normal during deformation. Whereas these assumptions 
simplify plate analysis, they also impose several restrictions on the form of 
polynomial approximations of the function w. In fact it can be shown that for 
monotonic convergence, it is necessary to match not only values of w at boundary 
points, but also slopes and twists at the corners of the element. With this in 
mind, the following polynomial approximation for w for a given element is 
introduced: 
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3 3 

w- 2 2 A (X,) r (X 2 ) 

r = 0 s = 0 


( 56 ) 


Here the 16 quantities A are undetermined constants. 


rs 

Let w N> e Na > and respectively denote the values of w, the first 

partial derivatives of w[9 = (w ) , T ] , and the mixed derivative of 

^ No: , a N 

w[£ = (w ) 1 at node N of a typical finite element. Then 

N , <2/3 N 


W N = TiTj A rs (X N1^ (X N2^ 

0 Ni = Z2 rA rs (X Nl ) (X N2 ) 
9 N2 = 22 sA rs (X Nl^ (X W 


and 


ZZ rsA 


(x M ,) r_1 <xi 

rs Nl N2 


s - 1 


Equations ( 57) represent 16 simultaneous equations in the 16 unknowns A . 

Solving these equations, introducing the results into equations (56) , and 
rearranging terms, 


w = H N W N + V ®N<* + J N hj 


(58) 


where the functions H XT , I XT , and J, T are defined in Table II. 

N No; N 

It should be noted that this displacement function is a two-dimensional 
generalization of the Hermite interpolation polynomial [1]. 
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TABLE H. POLYNOMIAL COEFFICIENTS - CASE I 


N 

1 

2 

3 

4 

N 

1 

2 

3 

4 

Note: 



Substituting equation { 58) into the strain-displacement relations (55) , 


2y ap 2Z(H N,Q!^ W N ^y.Qi/S 9 Ny + J N,Q!/3 


(59) 


The stiffness relation obtained by combining equations (59) and (39b) can 
be represented by a 16 x 16 bending-stiffness matrix. The bending-stiffness 
relation is of the form 




= k, u, 

~b Ho 


(60) 


where 

R, b = {Pl 3 » m il> m l2> qi;P23. m 2i • • • qj (6ia) 

U b = {w lf Qii> © 12 ) w 2 > ® 2 i> ••• ^4} ( 6 lb) 


and q N represents the twist applied at node N. This finite element formulation 

was first presented by Bogner, Fox, and Schmidt [11. The results obtained 
herein, which were derived independently, are in agreement with those found 
in Reference 1. 
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Because of the Kirchhoff hypothesis, the shear stiffness for this case has 
been assumed zero. 

Case n. The second representation studied assumes the transverse dis- 
placement w, and the rotations ©^ may be approximated by bilinear forms 


and 


0 

a 




(62) 


Note that this is a special case of equations (27) with i p = f^. 

The bending and shear stiffness relations obtained by substituting equa- 
tions (62) and (28) into equation (39) can be represented as 

p = (k, +k )U (63) 

~ ~b ~s ~ 


where, for this case, 

P = (Pi 3 » P 23 > P33 ’ P43; m ii> m 2l • • • m 42 ) 


(64a) 


and 


U = {w 1 , w 2 , Wj , w 4 ; ©u, © 21 , ... © 42 } • (64b) 

The elements of the bending and shear stiffness matrices obtained for 
this case are presented in Tables III and IV, respectively. 

Case III . It will be shown later that the bilinear approximation considered 
in Case II results in a representation that converges very slowly. To improve 
the convergence of the bilinear approximation, a discrete Kirchhoff hypothesis, 
which imposes constraints on the transverse shears, is introduced in this section. 
Specifically, it requires that the average transverse shear strain vanish along 
the edges. This is equivalent to requiring that the transverse shear vanish at 
the mid-point of each edge since the shear strains are linear along the edges. 

The geometry of deformation of a typical element edge subjected to this discrete 
Kirchhoff hypothesis is illustrated in Figure 3. 
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TABLE m. BENDING STIFFNESS MATRIX - CASE H 


1 2 
2(aj + a 2 ) (-2aj+a 2 ) 

2 ( aj + a 2 ) 


SYMMETRIC 


3 

4 


5 

6 


7 


-(aj + a 2 ) 

(aj - 

2 a 2 ) 

-( a 3 J - a 4 ) 

(a 3 - 

a -i) 

(a 3 + a 4 ) 

(a -l 

(a t - 2 a 2 ) 

-<»i 

a 2 ) 

( a 4 “ a 3) 

( a 3 ■ f * 

a 4> 

' a 3 - a 4 > 

-(a 

2 ( aj + a 2 ) 

(a 2 - 

2*0 

U 3 + a 4 ) 

(a 4 - 

a 3 ) 

-< a 3 - » 4 ) 

(a 3 


2 ( Qj 

+ a 2 ) 

( a 3 - a 4 ) 

-(a 3 

- a 4) 

(a 4 - ajl 

(a 3 




2 ( a j - (t 2 ) 

(°l - 

2d z ) 

-( O'l - o 2 ) 

l Cl 2 





2 (a - 1 

■** rt 2 > 

( Of 2 - 20!) 

-( ii 


2(0-4 i a 2 ) 


Note: Symbols used in this table are defined in the Appendix. 

TABLE IV. SHEAR STIFFNESS MATRIX - CASE II 

Row Column 

12345678 9 10 

1 16 3 4 8 0 0 0 0 12 d 5 6 d 5 

2 16 8 4 0 0 0 0 6 d 5 12 d 5 

3 16 8 0 0 0 0 6 d 5 12 d 5 

4 16 0 0 0 0 12 d 5 6 d 5 

5 16 8 4 8 -12 d, 12 d- 


16 8 


16 8 
16 


SYMMETRIC 


9 

10 

11 

12 d 5 

6 d 5 

-6 d 5 

6 d 5 

12 d 5 

-12 d 5 

6 d 5 

12 d 5 

-12 dj 

12 d 5 

6 d 5 

-6 d 5 

-12 d 4 

12 d 4 

6 d 4 

-12 d^ 

12 d 4 

6 d 4 

-6 d 4 

6 d 4 

12 di 

-6 d4 

6 

12 d 4 

18 d, 

9 d 2 

-9 d 4 


18 d 4 

9 d 3 



18 dj 


Symbols used in this table are defined in the Appendix. 


z 



(a) Undeformed State 


The discrete Kirchhoff hypo- 
thesis introduces a dependence between 

the rotations 0^ T and out-of-plane dis- 
Na 

placements w^. These relations may 
be expressed by 


Z 



U = 2b U* 

Aj 


(65) 


where U is defined by equation ( 64b) 
and U* is a vector of independent 
variables which contains eight compo- 
nents. An equivalent load vector g* 
may be obtained by equating the work 
done in the starred and unstarred sys- 
tem. 


FIGURE 3. GEOMETRY OF 
DEFORMATION OF A TYPICAL ELE- 
MENT EDGE UNDER THE DISCRETE 
KIRCHHOFF HYPOTHESIS 


*T * T . , 

p U = p U . (66) 

It follows from equations (64) and (65) 
that 

E* = & E • (67) 


A stiffness relation is obtained for the starred system by combining equations 
( 63) , ( 65) , and ( 67) . 


P* = (k* + k*) U* 

JQ g 


where 



and 

* T 

k* = ib k ip . 

oo g Ao ~ g Aj 


( 68 ) 


(69a) 


(69b) 


The elements of the bending and shear-stiffness matrices in the starred 
system are presented in Tables V and VI, respectively. In these tables the 
vector of independent variables is defined as 
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( 70 ) 


U* ={w 1 , w 2 w 3 , w 4 ; © llt 0 21 , Qje, © 42 } • 


TABLE V. BENDING STIFFNESS MATRIX - CASES III AND IV 


Row 

Column 


1 2 

3 

4 

5 

6 

7 

8 

1 

gi + 6A -gj 

Si - 6A 

-Si 

b S4 

bg 5 

-ag. 

-ag 5 


+2(g 2 +g 3 ) +g 2 - 2g 3 

‘§2 " S3 

-2g 2 + g 3 

+2/i a 1 a 

+2 \x a t a 

-4n« 2 b 

-2/i o? 2 b 

2 

gi - ex 

- Si 

Si 

-bg4 

-bg 5 

-ag 5 

-ag 4 


+2(g 2 + g 3 ) 

+ g 3 - 2g 2 

-<ga + 63 ) 

+ M 04 a 

+2/i c! t a 

+2/i cn 2 b 

+ /i of 2 b 

3 


gl + 6A 

“Si 

-bg 5 

-bg4 

*S 5 

a S4 



+ 2( g 2 + g 3 ) 

+ <g 2 - 2g 3 ) 

- p O'} a 

-2(jl a 2 a 

+ n a 2 b 

+2/i a 2 b 

4 



St - 6A 

b S 5 

bg 4 

ag 4 

a g5 




+ 2(g 2 + g 3 ) 

-2/i a 

- /x a? 1 a 

-p.oi 2 b 

-2/i a 2 b 

5 




Gb 2 

+4n a 2 

-Gb 2 
+2/i a 2 

-abA 

-abA 

6 





-Gb 2 
+4/i a 2 

-a b A 

-abA 

7 






Ga 2 
+4 fj. b 2 

-abA 

8 

1 







Ga 2 

+4 n b 2 

Note: 

Symbols used in this table are 

defined in the Appendix. 




Case IV . The main reason for including the transverse shear stiffness 
in the bilinear approximation is to allow treatment of problems involving forces 
acting normal to the middle surface. This results from the bending and mem- 
brance stiffness matrices being independent of the out-of-plane displacement w. 
Since the discrete Kirchhoff hypothesis, however, leads to a dependency between 
the rotations 9 and displacements w , it is no longer necessary to include the 

shear stiffness in the analysis. In addition, it is logical to omit the shear stiff- 
ness since the discrete Kirchhoff hypothesis causes the transverse shears to 
vanish in the limit. Therefore , in this case , the discrete Kirchhoff hypothesis 
is imposed on the bilinear approximation with the shear stiffness omitted. The 
stiffness relation is 



U* 


where k* is the same as in Case III. 
~b 


( 71 ) 
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TABLE VI. SHEAR STIFFNESS MATRIX - CASE HI 


Row 





5 

Column 

6 

7 

8 

1 

2 

3 

4 

1 

2 e 4 

e 2 


e 3 

2 e 6 

e 6 

-2 e 5 

-e 5 

2 


2 e 4 

e 3 

~ e i 

e 6 

2 e e 

2 e 5 

e 5 

3 



2 e t 

e 2 

~ e 6 

-2 e e 

e 5 

2 e 5 

4 




2 e 4 

-2 e e 

~ e 6 

-e 5 

-2 e 5 

5 





2 e 4 

e 4 

0 

0 

6 

SYMMETRIC 




2 e 4 

0 

0 

7 







2 e 4 

e 4 

8 








2 e 4 


Note: Symbols used in this table are defined in the Appendix. 


NUMERICAL EXAMPLES 
Square Plate Examples 


To test the merits of the stiffness matrices developed in this report, 
several numerical examples are presented in this Section. Four different plates 
are considered involving two edge support conditions and two loading conditions. 
The geometry of the plate and designations are shown in Figure 4. 

Because of the symmetry of this problem , only one quadrant of the plate 
was considered in the analysis. The convergence of each finite element re- 
presentation was obtained by considering different mesh sizes in the analysis of 
each case. Figure 5 illustrates the mesh arrangements studied. 
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Edge Length 
Thickness 



Edges: Simply Supported (SS) or Clamped (C) 

Loading: Uniformly Loaded (U) or Concentrated (C) 



Case 

Edges 

Load 

a 

C 

U 

b 

C 

C 

c 

SS 

U 

d 

SS 

C 




FIGURE 4. GEOMETRY OF PLATE FIGURE 5. TYPICAL FINITE ELE- 
AND CASE DESIGNATIONS MENT IDEALIZATION 

The numerical results obtained in this investigation are presented in 
Figures 6 through 9. The nondimensional central deflection is presented as a 
function of the number of elements in a quadrant of the plate. For the case of 
uniformly loaded plates, the dimensionless central deflection coefficient a is 
is given by 


a 


DW 
c 

qa 4 


(72) 


where D is the flexural rigidity, W is the central deflection, and q is the load 

c 

intensity. 


In the cases of the concentrated load, the deflection coefficient /3 is 
defined by 


J8 = 


DW 
c 

Pa 2 


(73) 
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where P is the concentrated load. The continuum analysis results are from 
Timoshenko and Woinowsky-Krieger [ li ] . Since Cases II and III include both 
bending stiffness, which is proportional to the thickness cubed, and shear 
stiffness, which is linearly proportional to the thickness, the results for these 
cases are dependent on the thickness-over-length ratio (h/a) . In this analysis 
the ratio h/a was set equal to 0. 1. 

The high-order polynomial approximation of Case I is seen to possess 
good convergence characteristics , as expected. These results , which were 
obtained independently, are in agreement with earlier results [1]. It should be 
emphasized that, because of the high-order polynomial approximation involved 
in this case, this representation proves cumbersome in the analysis of shells or 
geometrically nonlinear plate problems. A simple approximation is required 
to handle more complicated problems. 



FIGURE 6. CENTRAL DEFLECTION OF A SQUARE 
PLATE-CASE a ( C - U) 
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10 “ 



Number of Elements (NE) 

FIGURE 7. CENTRAL DEFLECTION OF A SQUARE PLATE-CASE b ( C - C) 



Number of Elements (NE) 


FIGURE 8. CENTRAL DEFLECTION OF A SQUARE PLATE-CASE c(SS-U) 




Number of Elements (NE) 

FIGURE 9. CENTRAL DEFLECTION OF A SQUARE PLATE-CASE d(SS-C) 

It is seen from the graphs that the simple bilinear approximation of Case 
II yields an element which is too stiff. This representation forces an unrealis- 
tically large portion of the total strain energy to be taken into shear energy. 
Usually in the analysis of thin plates the shear energy is negligible relative to 
the energy in bending. Therefore, the approximations which include the effects 
of shear, such as Cases II and III, result in an element which is too stiff. These 
approximations are more realistic for thicker plates. Consequently the numeri- 
cal results of Cases II and III will converge faster for thicker plates. 

The Case III representation results in an element which is even stiffer 
than that of Case II. In this approximation the shear is constrained to behave in 
a manner which is realistic for thin plates. Since, however, the shear energy 
is still included in the analysis, these constraints tend to make the element un- 
duly stiff. 

The Case IV approximation, which simply drops the shear stiffness from 
the Case III analysis, yields excellent results. These results are even more 
surprising because of the fact that this representation has only 16 degrees of 
freedom as compared with the 24 degrees of freedom of Case I. In addition, this 
representation uses a simple bilinear approximation. This finite element 
representation appears, therefore, to provide a simple basis for the treatment of 
plates and shells. 
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Frequently stresses must be computed in addition to the displacement. 
The stresses may be calculated from the nodal displacements by first applying 
the strain-displacement relations of equation ( 30) , and then the linear stress- 
strain relation of equation ( 48) . The finite element method described in this 
report is based on an assumed displacement pattern. Since the assumed dis- 
placement fields are continuous across element boundaries, include possible 
rigid body motions, and can lead to uniform strain states, the discrete model 
converges to the true continuum state of deformation as the network is refined. 
Stresses, however, converge in a mean square sense, but they do not, in 
general, converge monotonically. 

In this study the stresses were calculated at the center of a uniformly 
loaded simply-supported square plate. The stresses were computed at the four 
node points nearest the plate’s midpoint and then averaged. Figure 10 shows 
the nondimensional center stress (a/o^) and displacement (w/w q ) as a function 

of the number of elements, where a and w are the results obtained from a 
continuum analysis [11]. ° ° 



Number of Elements 


FIGURE 10. CENTRAL STRESS AND DEFLECTION OF A SQUARE 
SIMPLY-SUPPORTED UNIFORMLY LOADED PLATE 
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CONCLUSIONS 


The high-order polynomial approximation exhibits good convergence 
characteristics. The simple bilinear approximation with no constraints on 
the transverse shear strains yields an element that is too stiff. By introducing 
a discrete Kirchhoff hypothesis into the bilinear approximation and retaining both 
the bending stiffness and shear stiffness, an element is obtained which is even 
stiffer. When the shear stiffness is dropped from the discrete Kirchhoff bilinear 
representation, good convergence is obtained. This last approximation is the 
simple bilinear form and, in addition, contains only 16 degrees of freedom, 
whereas the first representation involves a high-order polynomial approximation 
with 24 degrees of freedom. Therefore, this finite element representation ap- 
pears to be readily applicable to the analysis of geometrically linear and non- 
linear plate and shell problems. 


George C. Marshall Space Flight Center 

National Aeronautical and Space Administration 

Marshall Space Flight Center, Alabama, June 28, 1968 
933-50-02-00-62 
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APPENDIX 


POLYNOMIAL COEFFICIENTS AND STIFFNESS MATRICES 


The polynomial coefficients and stiffness matrices developed in the text 
are presented in tabular form in this Appendix. Note that the symbols a, b, and 
h indicate the edge lengths in the X* and X 2 directions and the thickness of the 
plate, respectively. Young's modulus is denoted as E and Poisson's ratio, as 
v. A brief explanation of each table follows: 


Membrane Stiffness Matrix 


that 





~m 

_ J_ _b 



i - 3v 

6 a 


b 4 = 

8 

(1 -v) 

a 


1 a 

12 

b 

t>5 = 

6 b 

1 + v 


be = 

(1-v) 

4 n 


(A-i) 


and the matrix k is defined as 
~m 


T 

k « Kk 
~m ~m 


( A-2) 


where 


K = 


Eh 


( 1 - » 2 ) 


( A-3) 


and k is the matrix listed in Table I. 
~m 
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Polynomial Coefficients - Case I 


Table II lists the coefficients of the Hermitian interpolation polynomial 
used in the transverse displacements approximation of Case I. Recalling that 


w = H w + I 0 + J £ 

N N No Na N N 


(A- 4) 


the coefficients H^, I , and are listed for N = 1, 2, 3, 4, ce = 1, 2. For 
convenience the following nondimensional parameters are defined: 



and 


C 2 = 


X2 

b 


Bending Stiffness Matrix - Case 1 1 


The bending stiffness relation is given by 
&b = ~b ~b 


where 


u b = { 0 ». 


®21> ®31> ®41> ® 12 ’ ®22> ®32> 


and 


R b ={ m ii» m 21» m 31» m 41> m 12> m 22> m 32> m 42> • 

T 

The matrix listed in Table III is defined by the equation 

k = h3 k T 
~b 48 ab ~b ' 


(A-6) 


(A-7a) 


( A-7b) 


(A-8) 
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Table in uses the following nondimensional parameters: 


a l = 


Gb z 


a 2 = (X + 2G) — 


a 3 = X 


ab 


_ ab 
a 4 = G— 


(A-9) 


Q = 


Oi-> = — 

i a 


where 


G = 


E 


2(i + v) 


and 


X = 


Ev 


( 1 - v 2 ) 


( A-10) 


Shear Stiffness Matrix - Case 1 1 


The shear stiffness relationship is of the form 

p =k U , ( A-li) 

~s ~s ~s 

where 

Eg = m 2i> m 31> m 4i; m 12> m 22> m 32» Pl3> P23> P33’ P43} (A-12) 

and 
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u ={© 11 , ©21» ©32 > ©41’ °12> ©22 > ©32 ’ 9 42 > Wi, W 2 , W 3 , W 4 } 
s 


T 

The matrix k preseated for Case II in Table IV is defined as 
~s 


, Ghab , T 

k = — 777 k 
~s 144 ~s 


The following parameters are used in the table: 

_ 8(a 2 + b 2 ) 
dl “ 3 a 2 b 2 


_ 8( a 2 - 2b 2 ) 

2 3 a 2 b 2 

_ 8( b 2 - 2a 2 ) 

3 3 a 2 b 2 



(A-12b) 


(A-13) 


(A-14) 


Bending Stiffness Matrix - Cases 1 1 1 and IV 


The bending stiffness relationship for Cases in and IV is of the form 

p* = k* U* (A-15) 

K ~b ~ 


where 


U* = {Wi, W 2 ,W 3 , W 4 , ©u, © 2 i, ©a, © 42 ) . (A-16) 

and 

U = il) U* (A-17) 
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I 


where 


where 


and 


where 


U is defined by equation (A- 12b) . 
~s 

It follows that 


* ,T 

£ = £ Ec 


p is given by equation (A- 12a) . 

rog 

5}C 

The matrix presented in Table V is defined by 


k* _ u 

~b 36 ab ~b 


*T 


The following parameters are used in the table: 
gl = 14G 




g4 - 3A + G 
g 5 = 3X - G 



(A-18) 


(A-19) 


( A-20) 


H = A + 2G . 
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Shear Stiffness Matrix - Case 1 1 1 


The shear stiffness relationship for Case III is of the form 


p* = k* U* 

~ g 


( A-21) 


where U* and g* are defined by equations (A-16) and (A-18) , respectively. 
^ T 

The matrix k presented in Table VI is defined by 


k* 

~s 


Gh *T 
I8ab ~s 


( A-22) 


The following parameters are used in Table VI: 


e* = a 2 + b 2 

e 4 = a 2 b 2 


e 2 = a 2 - 2b 2 

e 5 = a b 2 

(A-23) 

e 3 = b 2 - 2a 2 

e 6 = a 2 b 
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